install libraries

library(raster)
## Warning: package 'raster' was built under R version 3.5.3
## Loading required package: sp
library(sp)
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.5.3
## rgdal: version: 1.4-2, (SVN revision 814)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/ASUS PC/Documents/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/ASUS PC/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
library(wesanderson)
## Warning: package 'wesanderson' was built under R version 3.5.3
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 3.5.2
library(tmap)
## Warning: package 'tmap' was built under R version 3.5.2

load data. raster must have background value of 0. clip by mask in arcmap.

b1 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF1.tif')
b2 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF2.tif')
b3 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF3.tif')
b4 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF4.tif')
b5 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF5.tif')
b6 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF6.tif')
b7 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF7.tif')

raster examination. sense check. vegetation should appear more bright in red and NIR.

par(mfrow = c(2,2))
plot(b2, main = "Blue", col = gray(0:100/100))
plot(b3, main = "Green", col = gray(0:100/100))
plot(b4, main = "Red", col = gray(0:100/100))
plot(b5, main = "NIR", col = gray(0:100/100))

stack false colour and rgb

landsatRGB <- stack(b3, b2, b1)
landsatFCC <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

## full landsat

files <- paste0('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF', 1:7, ".tif")
landsat_1981 <- stack(files)
names(landsat_1981) <- c('blue','green','red','nir','swir1','tir','swir2')
extent(landsat_1981)
## class       : Extent 
## xmin        : 252769.2 
## xmax        : 288859.2 
## ymin        : 522889.4 
## ymax        : 558259.4
landsat_81_sub <- subset(landsat_1981, 1:3)
nlayers(landsat_81_sub)
## [1] 3

band relationship

pairs(landsat_1981[[1:2]],main='blue v. green')

pairs(landsat_1981[[3:4]],main='red v. nir')

### normalised index

nd <- function(img, a, b){
  ba <- img[[a]]
  bb <- img[[b]]
  index <- (ba-bb)/(ba+bb)
  return(index)
}

### bare soil index. a= swir1, b= red, c= nir, d= blue
bli <- function(img,a,b,c,d){
  ba <- img[[a]]
  bb <- img[[b]]
  bc <- img[[c]]
  bd <- img[[d]]
  index <- ((ba+b3)-(bc+bd))/((ba+b3)+(bc+bd))
  }

ndvi, ndwi, ndbi

ndvi <- nd(landsat_1981, 4,3)
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')

ndwi <- nd(landsat_1981, 2,4)
plot(ndwi, col = rev(terrain.colors(10)), main = 'Landsat-NDWI')

ndbi <- nd(landsat_1981, 5,4)
plot(ndbi, col = rev(terrain.colors(10)), main = 'Landsat-NDBI')

ndbl <- bli(landsat_1981, 5,3,4,1)
plot(ndbl, col = rev(terrain.colors(10)), main = 'Landsat-NDBLI')

view histogram

hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab= "Frequency",
     col = "green",
     xlim = c(-1, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-1,1, 0.05), labels = seq(-1,1, 0.05))

hist(ndwi,
     main = "Distribution of NDWI values",
     xlab = "NDWI",
     ylab= "Frequency",
     col = 'steelblue3',
     xlim = c(-1, 0.8),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-1,0.8, 0.05), labels = seq(-1,0.8, 0.05))

hist(ndbi,
     main = "Distribution of NDBI values",
     xlab = "NDBI",
     ylab= "Frequency",
     col = "navajowhite",
     xlim = c(-1, 0.7),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-1,0.7, 0.05), labels = seq(-1,0.7, 0.05))

interquartile ranges

summary(ndvi)
##                 layer
## Min.    -4.857143e-01
## 1st Qu.  4.214876e-01
## Median   5.370370e-01
## 3rd Qu.  6.056338e-01
## Max.     7.815126e-01
## NA's     8.130520e+05
summary(ndwi)
##                 layer
## Min.    -6.666667e-01
## 1st Qu. -4.918033e-01
## Median  -4.368932e-01
## 3rd Qu. -3.543307e-01
## Max.     6.279070e-01
## NA's     8.130520e+05
summary(ndbi)
##                 layer
## Min.    -9.736842e-01
## 1st Qu. -2.456140e-01
## Median  -1.834320e-01
## 3rd Qu. -5.633803e-02
## Max.     5.902778e-01
## NA's     8.134540e+05

thresholding for land cover; cbind from min to first quartile

veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main = 'Veg cover')

water <- reclassify(ndwi, cbind(-Inf, -0.3, NA))
plot(water, main = 'Water cover')

built <- reclassify(ndbi, cbind(-Inf, -0.2, NA))
plot(built, main = 'Built cover')

sense check: overlay over RGB

plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Veg Over RGB")
plot(veg, add=TRUE, legend=FALSE)

plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Water Over RGB")
plot(water, add=TRUE, legend=FALSE)

plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Built Over RGB")
plot(built, add=TRUE, legend=FALSE)

export mask for arcmap classifier training

writeRaster(veg, filename= "1981_veg-mask.tif", overwrite=TRUE)
writeRaster(water, filename= "1981_water-mask.tif", overwrite=TRUE)
writeRaster(built, filename = "1981_built-mask.tif", overwrite=TRUE)